knitr::opts_chunk$set(echo = TRUE, digits = 3)
library(ggplot2)
library(dplyr)
سیستم تحلیل برآوردگرهای آماری
شبیه‌سازی مونت کارلو

1 مقدمه و چارچوب نظری

برآورد پارامترها یکی از چالش‌های اساسی در استنباط آماری است. در این پروژه، دو روش اصلی برآورد را مقایسه می‌کنیم:

  • روش گشتاوری (Method of Moments - MOM): گشتاورهای نمونه را با گشتاورهای جامعه برابر می‌گیرد
  • روش ماکسیمم درست‌نمایی (Maximum Likelihood Estimation - MLE): احتمال مشاهده داده‌ها را بیشینه می‌کند

هدف تحلیل:

در این تمرین، کارایی این دو روش را برای دو توزیع متفاوت (پارتو و پواسون) با استفاده از شبیه‌سازی مونت کارلو مقایسه می‌کنیم. معیارهای مقایسه شامل اریبی (Bias) و میانگین مربعات خطا (MSE) است.


2 بخش اول: توزیع پارتو

2.1 تعریف توزیع و پارامترها

توزیع پارتو با پارامتر شکل \(\alpha\) (مجهول) و پارامتر مقیاس ثابت \(x_m = 1\) دارای تابع چگالی احتمال زیر است:

\[ f(x; \alpha) = \begin{cases} \frac{\alpha}{x^{\alpha+1}} & x \geq 1 \\ 0 & x < 1 \end{cases} \]

تابع توزیع تجمعی (CDF):

\[ F(x) = 1 - \frac{1}{x^\alpha}, \quad x \geq 1 \]

کاربرد توزیع پارتو:

این توزیع برای مدل‌سازی داده‌های دم‌سنگین استفاده می‌شود، از جمله توزیع درآمد، اندازه شهرها، و قیمت سهام.

2.2 روش تبدیل معکوس برای شبیه‌سازی

برای تولید داده‌های تصادفی از توزیع پارتو، از روش تبدیل معکوس (Inverse Transform Sampling) استفاده می‌کنیم.

2.2.1 اشتقاق فرمول

با قرار دادن \(u = F(x)\) و حل معادله برای \(x\):

\[ u = 1 - \frac{1}{x^\alpha} \implies \frac{1}{x^\alpha} = 1 - u \]

\[ \implies x^\alpha = \frac{1}{1-u} \implies x = (1-u)^{-1/\alpha} \]

از آنجا که \(1-u\) نیز یکنواخت بین 0 و 1 است، می‌توان نوشت:

\[ x = u^{-1/\alpha} \]

2.2.2 تولید نمونه دستی

با فرض \(\alpha = 2\)، پنج عدد تصادفی یکنواخت تولید و به داده‌های پارتو تبدیل می‌کنیم:

شماره \(u_i\) محاسبه \(x_i = u_i^{-1/2}\) داده نهایی
1 0.64 \(0.64^{-0.5}\) 1.25
2 0.25 \(0.25^{-0.5}\) 2.00
3 0.16 \(0.16^{-0.5}\) 2.50
4 0.81 \(0.81^{-0.5}\) 1.11
5 0.04 \(0.04^{-0.5}\) 5.00

نمونه تصادفی شبیه‌سازی شده: \(\{1.25, 2.00, 2.50, 1.11, 5.00\}\)

2.3 استخراج برآوردگر گشتاوری (MOM)

برای توزیع پارتو با \(x_m = 1\)، امید ریاضی برابر است با:

\[ E[X] = \frac{\alpha}{\alpha - 1}, \quad \alpha > 1 \]

با برابر قرار دادن میانگین نمونه (\(\bar{X}\)) با میانگین نظری:

\[ \bar{X} = \frac{\alpha}{\alpha - 1} \]

حل برای \(\alpha\):

\[ \bar{X}(\alpha - 1) = \alpha \implies \bar{X}\alpha - \bar{X} = \alpha \]

\[ \implies \bar{X}\alpha - \alpha = \bar{X} \implies \alpha(\bar{X} - 1) = \bar{X} \]

\[ \implies \hat{\alpha}_{MOM} = \frac{\bar{X}}{\bar{X} - 1} \]

2.4 استخراج برآوردگر ماکسیمم درست‌نمایی (MLE)

تابع درست‌نمایی برای نمونه \(x_1, \ldots, x_n\):

\[ L(\alpha) = \prod_{i=1}^{n} \frac{\alpha}{x_i^{\alpha+1}} = \alpha^n \prod_{i=1}^{n} x_i^{-(\alpha+1)} \]

لگاریتم درست‌نمایی:

\[ \ln L = n\ln(\alpha) - (\alpha+1)\sum_{i=1}^{n}\ln(x_i) \]

مشتق نسبت به \(\alpha\) و برابر صفر قرار دادن:

\[ \frac{d\ln L}{d\alpha} = \frac{n}{\alpha} - \sum_{i=1}^{n}\ln(x_i) = 0 \]

\[ \implies \hat{\alpha}_{MLE} = \frac{n}{\sum_{i=1}^{n}\ln(x_i)} \]

مقایسه فرمول‌ها:

  • MOM: \(\hat{\alpha}_{MOM} = \frac{\bar{X}}{\bar{X} - 1}\) — فقط به میانگین وابسته است
  • MLE: \(\hat{\alpha}_{MLE} = \frac{n}{\sum\ln(x_i)}\) — از اطلاعات لگاریتمی استفاده می‌کند

2.5 شبیه‌سازی مونت کارلو برای پارتو

set.seed(123) n <- 20 alpha <- 2 xm <- 1

set.seed(123)
n <- 20
alpha <- 2
xm <- 1

mom <- numeric(100)
mle <- numeric(100)

for (i in 1:100) {
  u <- runif(n)
  x <- xm / (u^(1/alpha))
  
  x_bar <- mean(x)
  mom[i] <- x_bar / (x_bar - xm)
  
  log_sum <- sum(log(x / xm))
  mle[i] <- n / log_sum
}

2.6 1. محاسبه اریبی (Bias)

 mean(mom) - alpha
## [1] 0.1658767
 mean(mle) - alpha
## [1] 0.05332796

2.7 2. محاسبه میانگین مربعات خطا (MSE)

 mean((mom - alpha)^2)
## [1] 0.1906337
 mean((mle - alpha)^2)
## [1] 0.1349972

2.8 نمودار مقایسه‌ای پارتو

رسم نمودار جعبه‌ای برای مقایسه چشمی

boxplot(mom, mle, names=c("MOM", "MLE"), col = c("lightblue", "lightgreen"),
main="مقایسه برآوردگرها برای توزیع پارتو",
ylab="مقدار برآورد α")
abline(h = alpha, col = "red", lty = 2, lwd = 2)
legend("topright", legend="α واقعی = 2", col="red", lty=2, lwd=2)

تحلیل نتایج پارتو:

  • اریبی: MLE اریبی کمتری دارد (~0.069) در مقایسه با MOM (~0.164)
  • MSE: MLE با MSE کمتر نسبت به MOM دقیق‌تر است
  • پراکندگی: هر دو روش واریانس مشابهی دارند اما مرکز توزیع MLE به مقدار واقعی نزدیک‌تر است
  • نتیجه: برای توزیع پارتو، MLE روش برتر است

3 بخش دوم: توزیع پواسون

3.1 تعریف و ویژگی‌های توزیع

متغیر تصادفی \(X\) با توزیع پواسون (نرخ \(\lambda\)) دارای تابع احتمال زیر است:

\[ P(X=x) = \frac{e^{-\lambda} \lambda^x}{x!}, \quad x = 0, 1, 2, \ldots \]

ویژگی کلیدی: پارامتر \(\lambda\) هم میانگین و هم واریانس توزیع است:

\[ E[X] = \text{Var}(X) = \lambda \]

کاربرد توزیع پواسون:

مدل‌سازی تعداد رخدادها در بازه زمانی ثابت: ورود مشتریان، تماس‌های تلفنی، خرابی دستگاه، و غیره.

3.2 استخراج برآوردگر گشتاوری (MOM)

از آنجا که \(E[X] = \lambda\)، برآوردگر گشتاوری به سادگی میانگین نمونه است:

\[ \hat{\lambda}_{MOM} = \bar{X} = \frac{1}{n}\sum_{i=1}^{n}x_i \]

3.3 استخراج برآوردگر ماکسیمم درست‌نمایی (MLE)

تابع درست‌نمایی:

\[ L(\lambda) = \prod_{i=1}^{n} \frac{e^{-\lambda} \lambda^{x_i}}{x_i!} = \frac{e^{-n\lambda} \lambda^{\sum x_i}}{\prod x_i!} \]

لگاریتم درست‌نمایی:

\[ \ln L = -n\lambda + \left(\sum x_i\right) \ln(\lambda) - \ln\left(\prod x_i!\right) \]

مشتق نسبت به \(\lambda\):

\[ \frac{d\ln L}{d\lambda} = -n + \frac{\sum x_i}{\lambda} = 0 \]

\[ \implies \hat{\lambda}_{MLE} = \frac{\sum x_i}{n} = \bar{X} \]

نتیجه مهم:

برای توزیع پواسون (و همچنین توزیع نرمال و دوجمله‌ای)، روش گشتاوری و ماکسیمم درست‌نمایی به فرمول یکسانی می‌رسند:

\[\hat{\lambda}_{MOM} = \hat{\lambda}_{MLE} = \bar{X}\]

این یعنی در شبیه‌سازی، عملکرد دو روش کاملاً یکسان خواهد بود.

3.4 مثال محاسباتی دستی

با \(\lambda = 4\) و نمونه \(\{3, 5, 2, 6, 4\}\):

\[ \bar{X} = \frac{3+5+2+6+4}{5} = \frac{20}{5} = 4 \]

بنابراین: \[ \hat{\lambda}_{MOM} = \hat{\lambda}_{MLE} = 4 \]

این برآورد دقیقاً با مقدار واقعی پارامتر برابر است (که البته تصادفی است).

3.5 شبیه‌سازی مونت کارلو برای پواسون

n <- 30
lambda <- 4

mom <- numeric(100)
mle <- numeric(100)

for(i in 1:100){
  x <- rpois(n, lambda)
  
  mom[i] <- mean(x)
  mle[i] <- mean(x) 
}
mean(mom) - lambda
## [1] -0.02633333
mean(mle) - lambda
## [1] -0.02633333
mean((mom - lambda)^2)
## [1] 0.1317667
mean((mle - lambda)^2)
## [1] 0.1317667

3.6 نمودار مقایسه‌ای پواسون

رسم نمودار جعبه‌ای

boxplot(mom, mle, names=c("MOM", "MLE"), col = c("orange", "orange"),
        main="مقایسه برآوردگرها برای توزیع پواسون",
        ylab="مقدار برآورد λ")

abline(h = lambda, col = "red", lty = 2, lwd = 2)

legend("topright", legend="λ واقعی = 4", col="red", lty=2, lwd=2)

# خط زیر اصلاح شد: افزودن تابع text و بستن پرانتز
text(x = 1.5, y = lambda + 0.4, labels = "توجه: دو جعبه کاملاً منطبق هستند", 
     col = "darkblue", cex = 0.8)

تحلیل نتایج پواسون:

  • همسانی کامل: هر دو روش عملکرد یکسانی دارند (اریبی و MSE برابر)
  • اریبی ناچیز: اریبی هر دو روش تقریباً صفر است
  • MSE پایین: MSE نشان‌دهنده دقت بالا با n=30 است
  • نتیجه: برای پواسون، انتخاب روش اهمیتی ندارد

4 نتیجه‌گیری نهایی

4.1 مقایسه جامع دو توزیع

خلاصه دستاوردها:

  1. توزیع پارتو: MLE برتری قابل توجهی دارد با اریبی و MSE کمتر. دلیل: استفاده بهینه از اطلاعات لگاریتمی و مقاومت بهتر در برابر دم‌سنگینی
  2. توزیع پواسون: هیچ تفاوتی بین دو روش وجود ندارد زیرا هر دو به فرمول یکسانی (میانگین نمونه) می‌رسند
  3. انتخاب روش: برای توزیع‌های پیچیده (مانند پارتو)، MLE معمولاً کارآمدتر است، اگرچه محاسباتی سنگین‌تر است
  4. تأیید تئوری: نتایج شبیه‌سازی مونت کارلو کاملاً با پیش‌بینی‌های تئوری مطابقت دارد

کاربردهای عملی:

  • مالی: تحلیل داده‌های مالی با توزیع پارتو (درآمدها، ریسک)
  • صنعتی: تحلیل صف و تعداد خرابی‌ها (پواسون)
  • بیمه: مدل‌سازی فراوانی خسارت (پواسون) و شدت خسارت (پارتو)
  • رگرسیون: برآورد پارامترها در مدل‌های خطی تعمیم‌یافته

4.2 تحلیل نظری

چرا MLE برای پارتو بهتر است؟

  1. استفاده کامل از اطلاعات: MLE از تابع درست‌نمایی کامل استفاده می‌کند، در حالی که MOM فقط از گشتاور اول بهره می‌برد
  2. مقاومت به دم‌سنگینی: در توزیع پارتو، میانگین حساس به مقادیر افراطی است، اما MLE با لگاریتم این حساسیت را کاهش می‌دهد
  3. کارایی مجانبی: طبق قضیه کرامر-رائو، MLE در نمونه‌های بزرگ به کارآمدترین برآوردگر میل می‌کند

چرا در پواسون تفاوتی نیست؟

توزیع پواسون از خانواده توزیع‌های نمایی است و دارای آماره کافی (Sufficient Statistic) ساده‌ای به نام \(\sum X_i\) است. در چنین مواردی، هر روش معقول برآوردی به همان آماره کافی می‌رسد و عملکرد یکسانی دارند.